#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include <time.h>
#include <getopt.h>
#include <omp.h>

static int SIZE;
static int NUM_THREAD;

double **gMatrixA;
double **gMatrixB;

struct option input_option[] = {
	{"size", 1, 0, 's'},
	{"thread_count", 1, 0, 't'},
	{0, 0, 0, 0}
};

struct timeval startTime, endTime;

#if 0 
// For testing
void showMatrix()
{
	int i, j;
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE; ++j)
			printf("%f\t", gMatrixA[i][j]);
		printf("\n");
	}
	printf("\n");
}
#endif

int init()
{
	int i, j;

	/* Memory allocation */
	gMatrixA = (double**)malloc(sizeof(double *)*SIZE);
	for(i=0; i<SIZE; ++i)
		gMatrixA[i] = (double*)malloc(sizeof(double)*SIZE+1);

	gMatrixB = (double**)malloc(sizeof(double *)*SIZE);
	for(i=0; i<SIZE; ++i)
		gMatrixB[i] = (double*)malloc(sizeof(double)*SIZE+1);

	srand48(time(0));
	for(i=0; i<SIZE; ++i) {
		for(j=0;j<SIZE+1; ++j)
			gMatrixA[i][j] = gMatrixB[i][j] = drand48();
	}
	return 0;
}

int fini()
{
	int i;
	for(i=0; i<SIZE; ++i) {
		free(gMatrixA[i]);
		free(gMatrixB[i]);
	}
	free(gMatrixA);
	free(gMatrixB);

	return 0;
}

int gaussian_elim()
{
	int i, j, k;
	double dCoeffi;
	int chunk_size = SIZE / NUM_THREAD;

	#pragma omp parallel private(i, j, k, dCoeffi) shared(gMatrixA, chunk_size) num_threads(NUM_THREAD)
	{
		#pragma omp master
		gettimeofday(&startTime, NULL);

		for(i=0; i<SIZE; ++i){
			#pragma omp for schedule(static, chunk_size)
			for(j=0; j<SIZE; ++j){
				if(i==j)	continue;

				dCoeffi = gMatrixA[j][i]/gMatrixA[i][i];
				for(k=0; k<SIZE+1; ++k){
					gMatrixA[j][k] = gMatrixA[j][k] - dCoeffi*gMatrixA[i][k];
				}
			}
		}
	}

	#pragma omp parallel for schedule(static, chunk_size) private(j) shared(gMatrixA, chunk_size) num_threads(NUM_THREAD)
	for(j=0; j<SIZE; ++j){
		gMatrixA[j][SIZE] /= gMatrixA[j][j];
		gMatrixA[j][j] = 1.0;
	}

	return 0;
}

int getL2Norm(double* pL2Norm)
{
	double dL2Norm=0;
	double dTmp;
	int i, j;

	for(i=0; i<SIZE; ++i){
		dTmp = 0;
		for(j=0; j<SIZE; ++j){
			dTmp += gMatrixB[i][j]*gMatrixA[j][SIZE];
		}
		dTmp -= gMatrixB[i][SIZE];
		dL2Norm += pow(dTmp, 2.0);
	}
	*pL2Norm = sqrt(dL2Norm);

	return 0;
}

void set_option(int argc, char** argv)
{
	int opt;
	if(argc!=3 && argc!=5 ){
		printf("Usage: %s [options]\n", argv[0]);
		printf("\t-s\t--size\t\t Size of Matrix\n");
		printf("\t-t\t--thread_count\t Number of Thread\n");
		exit(0);
	}

	do{
		opt = getopt_long(argc, argv,"s:t:", input_option, NULL);

		switch(opt)
		{
			case 's':
				SIZE = atoi(optarg);
				break;
			case 't':
				NUM_THREAD = atoi(optarg);
				break;
			case -1:
				break;
			default:
				abort();
		}
	}while(opt!=-1);
}

int main(int argc, char** argv)
{
	double dL2Norm;

	set_option(argc, argv);

	init();
	
	gaussian_elim();
	gettimeofday(&endTime, NULL);

	getL2Norm(&dL2Norm);
	printf("\nL2Norm: %f\n", dL2Norm);
	printf("Time Ellapsed: %f msec\n", (endTime.tv_sec - startTime.tv_sec)*1000 + (endTime.tv_usec - startTime.tv_usec)/1000.0);

	fini();

	return 0;
}
